function [a,flag,param,N] = BoltzmannE(param, Q, aFeasible, iState, aFCat, a, t, n, N, MonthNumber)
% implementation of Boltzmann exploration
% see Powell, Approximate Dynamic Programming (2007), pp. 328-329
% niels riegels, 01.10.08

%compute beta parameter
if nargin<10, MonthNumber = 1; end

te = Q(iState,1:aFCat,MonthNumber); 
iAMax = find(te == max(te));
iAMax = iAMax(1); % get rid of non-uniqueness
iAMin = find(te == min(te)); 
iAMin = iAMin(1); % get rid of non-uniqueness
dQ = Q(iState,iAMax,MonthNumber) - Q(iState,iAMin,MonthNumber);
if dQ == 0
    dQ = 1e+3;
end
N(iState,1) = N(iState,1)+1;
beta = N(iState)/dQ;

%compute probability density function and cumulative distribution
P=beta*te;
P=exp(P);
totalP=sum(P);
if totalP
    P=P./totalP;
    C = cumsum(P);
end

%select action
rN=rand;
if ~totalP
    action=ceil(rN*aFCat);
else
    action=sum(C<=rN)+1;
end
a(t) = param.actionV{n}(action);
flag=0;